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Fundamental properties of superfluids with d-wave pairing symmetry are investigated theo- 
retically. We consider neutral atomic Fermi gases in a harmonic trap, the Cooper pairing being 
produced by a Feshbach resonance via a d-wave interaction channel. A Ginzburg-Landau (GL) 
functional is constructed which is symmetry constrained for five component order parameters 
(OPs). We find stable OP textures and vortices for all the three phases which are known to 
be the energy minimum of the GL functional; the ferromagnetic, polar and cyclic phases both 
at rest and under rotation. In particular, we touch upon the stability conditions for a non- 
Abelian fractional 1/3- vortex in the cyclic phase under rotation. It is proposed how to create 
the intriguing 1/3- vortex experimentally in atomic gases via optical means. 
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1. Introduction 

Superfluids with s-wave pairing symmetry have been 
realized by using a Feshbach resonance of ^Li atom gases 
in 2005^) by sweeping a magnetic field. It is natural to ex- 
pect that the research front of cold atom gases develops 
towards reafizing condensates with a higher relative an- 
gular momentum of the Cooper pairs. In fact, much effort 
from theoretical and experimental sides is now focused 
on p-wave pairing in ^Li where a p-wave Feshbach reso- 
nance is already confirmed to exist. Although p-wave 
molecules have been formed, experimental evidence of 
p-wave superfluidity in ^Li has not been reported yet. 

It might be useful to theoretically investigate super- 
fluidity with higher angular momentum / of the Cooper 
pair, such as d-wave {I = 2), or /-wave (/ = 3), in neu- 
tral Fermi gases to stimulate experimental and theoreti- 
cal works towards this direction. Hulet has performed a 
coupled channels calculation and found a Feshbach res- 
onance for d-wave channel in ^Li,'^) and hence we have 
a good reason to explore this possibility from a theoret- 
ical point of view. There have already been appearing 
theoretical studies towards this direction.^' 

It is known that the pairing symmetry of high Tc 
cuprates is described by the dx^-y^ state. The study of 
superfluids with d-wave symmetry might be important 
because several strongly correlated superconductors, or 
so-called heavy Fermion superconductors belong to un- 
conventional pairing states, such as d, d -I- id, or /, etc, 
including high T^ whose pairing mechanism is still un- 
clear. Note that in the heavy Fermion superconductor 
UPta a /-wave pairing is realized.^"*-' It is also inter- 
esting because of the richer physics associated with the 
many internal degrees of freedom of the relative orbital 
angular momentum / > 1 of the Cooper pairs compared 
to the / = (s) and I = 1 {p) cases. In particular, the 
spatial structure of the order parameter (OP), or the tex- 
tures and vortices, which are hallmarks of superfluidity, 
are intriguing in the confined ultracold atomic gases. 

We note that the OP structures of the d-wave su- 
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perfluid of Fcrmionic atoms and spinor condensates of 
Bosonic spin-2 atoms are mathematically similar because 
both OPs have 5 components. Whereas the latter has 
been extensively investigated already, ^^^^^ the former has 
not been studied so far in connection with neutral atom 
gases. The boundary conditions due to the harmonic trap 
are different between the two cases in an essential way, 
giving rise to novel textures and vortices in the d-wave 
superfluid. This difference arises because of the different 
origin of the OP degeneracy due to the internal degrees 
of freedom: The orbital angular momentum lives in real 
space in the d-wave superfluid unlike the spin degrees 
of freedom in the Bosonic case although both are of the 
same 5*0(3) symmetry in a homogeneous system. 

Mermin^^^ presents a general framework based on 
Ginzburg-Landau (GL) formalism to describe a d-wave 
superfluid in an infinite bulk system. He exhaustively 
classifies the ground state phase diagram into three re- 
gions; ferromagnetic (FM), polar (PO) and cyclic (CY) 
phases, depending on the coupling constants or fourth or- 
der coefficients /3i, (32, and /33 in the GL functional. The 
resulting phase diagram in the (/3i, /Js) plane normalized 
by (32 is shown in Fig. 1 where the three phases are enu- 
merated: the ferromagnetic (FM), polar (PO), and cyclic 
(CY) phases. The cyclic phase is the most intriguing one 
among them because Fermion superfiuids with s-wave 
(/ = 0) or spinless p-wave {1 = 1) pairings produced by a 
Feshbach resonance via a magnetic sweep, or the spin-1 
spinor BEC^''' do not support this phase which only 
appears above I > 2. 

The motivation of this paper is to provide the funda- 
mental physical properties of a d-wave superfluid con- 
fined by a harmonic potential in two-dimensions (2D) 
in order to help to identify the d-wave nature experi- 
mentally, focusing on the OP spatial structures or tex- 
tures at rest and vortices under rotation. We study the 
most stable textures and vortices for each phase under 
trapped geometry. In particular, the non-Abelian 1/3- 
vortex which is the energetically favored state under ex- 
ternal rotation in the cyclic phase is examined in some 
details. The existence of a 1/3-vortex provides a spectro- 
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Fig. 1. (Color online) Phase diagram of the stable phases in uni- 
form system where /3i vs /33 plane normalized by /32 is shown; 
ferromagnetic (FM), polar (PO) and cyclic (CY) phases, (a) cor- 
responds to the weak coupling case, (b), (c) and (d) correspond 
to the parameter values adopted in the strong coupling cases. 

scopic means to characterize d-wave superfluidity. Mu- 
tually non-commutative 1/3- vortices themselves are al- 
ready pointed out to be allowed topologically in an = 2 
spinor BEC by several authors. However, there is no 
serious calculation to examine its stability from energetic 
point of view in d-wave Fermion superfluids, which is one 
of our main purposes in this paper. As for the n F — 2 
spinor BEC we studied the energetics of 1/3- vortices in 
our recent paper. The non-Abelian 1/3-vortices might 
turn out to be useful for quantum computation because 
a spatial arrangement of mutually non-commutative 1 /3- 
vortices can be used for storing information, namely they 
can be used as a topologically protected qubit. This is 
somewhat similar to the idea that utilizes the Majorana 
particles bound to the vortex core in chiral p + ip super- 
conductors.^"^ Note that the 1/3- vortex discussed here 
does not support Majorana zero energy particles. 

The arrangement of this paper is following. In §2, we 
give a detailed derivation of the GL functional form rel- 
evant to the present trapped systems. The GL gradient 
term, which is important in the trapped Fermi gases, 
is derived here. We examine the case in which the GL 
coefficients are given by weak coupling estimate. This 
case is situated at the boundary between FM and CY in 
the fourth order parameter space mentioned above (see 
(a) in Fig. 1), finding that the FM phase is stabilized in 
trapped systems over the CY phase. In §3 we examine the 
spatial structures, or textures of the three phases when 
trapped, and also vortices under rotation for each phase 



within the strong coupling case. In particular, the 1/3- 
vortex in CY is discussed in detail because it is unique in 
CY phase and they are topologically interesting due to 
their non-commutative properties. In the final section we 
will devote to conclusions and summary and touch on an 
experimental method to create a 1/3-vortex in neutral 
Fermi atomic gases. A part of this paper is published in 
ref. 21. 

2. Formulation and Preliminary 

2. 1 Ginzburg- Landau functional 

Our arguments are based on GL theory which is gen- 
eral enough to describe a d-wave superfiuid in terms of 
the free energy expanded up to fourth order in the OP. 
GL theory is valid near the transition temperature Tc, 
however, its range of applicability is known to be wider 
empirically than that. The OP A(r, k) for a d-wave su- 
perfiuid is spanned by the spherical harmonics Yim (k) (k 
is a unit vector in the momentum space.) 

A(r,k) = Ara{v)Yi^{i^) = k,B,,k, (1) 

with ; = 2, m = -2, • • • ,2 and i,j = 1,2,3. The re- 
peated indices are summed over. The coefficient Am(r) 
is a complex valued functions of r = {x,y). To empha- 
size the essential points, we consider a two-dimensional 
system assuming that the obtained objects extend uni- 
formly to the third dimension. Alternatively to Am we 
may sometimes use the symmetric traceless 3x3 matrix 
Bij . The i?- matrix that has five independent elements is 
convenient in constructing the GL functional for the OP 
with 5*0(3) symmetry in addition to U{1) gauge sym- 
metry. In the following we use these two notations inter- 
changeably. The bulk GL free energy functional fbuik is 
derived by Mermin^'^^ as 

fbuik^aTrB*B + Pi\TrB^\'' 

+ P2{TrB*Bf + p^TriB^B*^) (2) 

where a{T) = ao{T — Tc) and Tc is the transition tem- 
perature. There are three independent fourth order terms 
/?!, P2 and /Sa. The trace operation for a matrix is denoted 
by Tr. After some algebra this bulk energy fbuik is recast 
into 

hulk = ao(T- T,)S,|A,|2 + y |/32 + 

+ (/3i + ^Psm' - ^/33|/T} V.f (3) 

where the orbital singlet pairing amplitude 

e^i-iyA,A^,/^i\Ai\^ (4) 

and the orbital momentum 

f^A*F,,A,/^i\Ai\' (5) 

with F being the 5x5 spin matrix. 

The three phases are characterized by FM (|0| = 
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0, I/I 0), CY (|e| - 0, 1/1 - 0), PO ([91 ^ 0, 1/1 = 0). 
It is clear from eq. (3) that CY is stable when /3i+/33/6 > 
and /33 < 0, while FM and PO phases occupy the other 
regions of the (/33,/3i) parameter space normalized by 
/32 (see Fig. 1). The canonical CY phase is of the form 
A(k) = iy22(k) + V2Y2o(k) + iY2-2^) or in a vecto- 
rial notation: (A2, Ai, Aq, A_2)^ = (i, 0, ^/2, 0, z)"^. 
The other cyclic states are obtained from this by rota- 
tions. The weak coupling estimate for the Fermi sphere 
(/32 = 2/3i and = 0) predicts that the stable phase is 
on the boundary between the CY and FM phases (see 
(a) in Fig. 1). 

The gradient terms are constructed by considering the 
possible contractions of (diBjk)* and (diBmn)] namely 
there are three independent terms: (1) (diBjk)* (diBjk) 
(2) (diBjk)* (dkB.j) and (3) (a,Bi,-)*(5feAfc)- This can be 
understood from the fact that the decomposition of the 
two angular momenta (L = 1) x (i = 2) ^ L = 3, 2, 1 
where the former (latter) corresponds to d (OP). Then, 
the gradient energy^^' ^"'^ can be summed up as 

fgrad = Ki{dtB jk)* {OiB jk) + K2{diBjk)* {dkB.j) 

+ Ks{djBijndkB,k) (6) 

where Ki = K, K2 = 2K, K3 = 2K with K = 2Ko/35 
in the weak coupling approximation which also gives 
a = 200/15^/32 = 2/3i = 8/?o/315, OP a mplitude 
A() = ^/ao/2f3o and the coherence length ^0 — v^o/oo- 
For ao = Np/Tc, /3o = 7C{3)NF/16Tr'^klT^ , and Kq = 
7C{S)NF{hvFf /'iS-rr^klT^. We can derive the result of 
eq. (6) also starting from the A^, representation after 
simple but tedious calculations, which is described in Ap- 
pendix in detail. 

In addition to these bulk and gradient terms, we take 
into account a harmonic trap potential: 

f harmonic — CJ T Xjj|j4.^| . (7) 

External rotation may be included by replacing d ^ d — 
i{2ni/h)fl X r where the rotation axis fl || z. Therefore 
the resulting total GL functional is given by 

f total — fbulk ~t- fgrad ~t- f harmonic + fceni (8) 

where fccni. is the centrifugal potential. The associ- 
ated mass current is derived as 

Arn 

ji = —lm[KrB*^ViBjk + K2B*^VkBij 

+ K:,B*^\7,B,,,]. (9) 

The following numerical calculations are done in two- 
dimensional plane of Cartesian coordinates {x,y). The 
coiiplcd GL equations for five components are solved on 
discretized lattice using i6ixl61 meshes via a simple it- 
eration. The OP amplitudes and length are scaled by Ao 
and ^0 respectively. The trap frequency uj is normalized 
by I '2-^^ ■^it;h m the atomic mass. We fix w = 1/80 in 
these units throughout this paper (except for the cases 



explicitly stated otherwise). The temperature is fixed to 
T = 0.8Tc. The field of view in the OP contour maps is 
30^0 X 30^0- The rotation frequency fl is normalized by 
the trap frequency ui. 

2.2 Weak coupling calculation 

Here before discussing the detailed results on each 
phase in the strong coupling case, we briefiy touch on 
the weak coupling case, which indicates the stable phase 
on the boundary between FM and CY for an infinite 
system as shown in (a) of Fig. 1, in order to understand 
the importance of the boundary condition through the 
gradient terms. We perform numerical calculations for a 
trapped system by solving the GL equations for five com- 
ponents. It turns out that FM is energetically favorable 
over CY in the weak coupling case. This is because of 
the boundary condition, which strongly constraints the 
possible phase. In fact, FM is more flexible than CY in 
the sense that the /-vectors can arrange to stabilize the 
state. As will be seen later (see Fig. 3 (a)), the /-vectors 
tend to align along the circumference on the boundary. 
This arrangement is energetically advantageous because 
the four point nodes situated along the / direction, which 
is explained later in §3.1, tend to lie outside the system, 
leading to maximal gain in the condensation energy. This 
is analogous to that in p-wave 

superfluids.25.26) This flex- 
ible feature is absent in CY because |/| = 0. Thus un- 
der confinement the degeneracy between FM and CY is 
lifted, the latter being never stabilized over FM in the 
weak-coupling limit. 

3. Strong coupling calculation 

3.1 Ferromagnetic phase 

We start out with the FM phase by choosing the pa- 
rameter values /3i = 1 and Pa — I (see (b) in Fig. 1). 
In Fig. 2 we show the contour maps of the OP magni- 
tude, I/I, and |e| for throe cases: (a) n = 0, (b) fl = 0.2, 
and (c) f2 = 0.25. It is seen from those figures that at 
rest the moment vector lies in the 2D plane, pointing to 
the horizontal .r-direction mostly and curving along the 
circumference regions. There is no fz component at rest 
as seen Fig. 3(a). The non- vanishing components consist 
exclusively of Ai, Aq, and A^i as seen from Fig. 2(a), 
which form the in-plane moment components fx and fy. 
Numerical data show that in our solutions the relation 

Aj = A*_j (10) 
holds at every spatial point at rest. Since 

/, = (2IA2P + |Ai|2 - |A_i|2 - 2|A_2n/|A|2 (11) 
this leads to /2 = 0. As for 6 it is now given by 

e ={2\A2\^-2\Ai\^+Al)/\A\\ (12) 

Because of the fact that |^2|^ = |^-2p < 1, 6 = 
{Al - 2|^i|2)/|A|2. Thus if |Ao| = V2\Ail 9 = 0. Ac- 
cording to our numerical data the above is nearly satis- 
fied. Therefore, the obtained phase shown in Fig. 2(a) is 
genuine ferromagnetic state. 
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Fig. 2. (Color online) Contour maps of IA2I ~ |^-2|, |/|, and |e| in FM phase. /3i = 1 and /Ss = 1. (a) = 0, (b) f2 = 0.2, and (c) 
n = 0.25. 
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Fig. 3. (Color online) The fx-fy vectorial patterns are shown, 
corresponding to (a), (b) and (c) respectively in Fig. 2. The arrow 
(color code) denotes fx-fy (fz) components. 



As for the nodal structure of this FM phase, we start 
with the most general expression for the c?-wave OP: 



A2V3{k + iky? - Ai2V3{k 
-Ao\/2(3fc2 - 1) + A_i2V3(A:^ 



iky)kz 
(13) 



The nodal condition |A| — is generally satisfied by 
the four point nodes on the Fermi sphere. At the central 
region where / aligns horizontally the four point nodes 
situated at 6* - 90° ± 15° and (^=0° and 180° as shown 
in Fig. 4(a) where the Fermi sphere is described by 9 
and (f) in the spherical coordinates. At the circumference 
region the /-vector tends to align so that two of the four 
point nodes are outside of the condensate. We display 
an example in Fig 4(b) where the two point nodes at 
9 = 90° and (f) '^130° and 315° are nearly outside the 
condensate. Therefore, it is understood that the curving 
of the /-vectors in the boundary region is due to the 
saving of the condensation energy by excluding the point 
nodes from the condensate. This situation is the same as 
in the p-wave condensate in the trapped systems. ^^'^^^ 

Under rotation Vl = 0.2 shown in Fig. 2(b) and 
Fig. 3(b) the A^i component exhibits phase singulari- 
ties on the y-axis seen as two dots in Fig. 2(b). Those 




Fig. 4. (Color online) Order parameter amplitude |A| in fc-space 
represented by the spherical coordinates {9, <j)) on the unit Fermi 
sphere, (a) At the central position and (b) at the upper right 
corner in Fig. 3(a). The darker regions indicate the nodal portion 
of 1A|. 



singularities in A_i are filled by the Ai component seen 
as two yellow points. Likewise the two singularities in the 
A-2 component are filled by the A2 component. At those 
singularities the /-vectors wind around near those cores 
and fz component is induced as seen from Fig. 3(b). In 
addition, the other components such as A±2 are induced 
by rotation (see Fig. 2(b)). Upon further increasing the 
rotation to il = 0.25, Fig. 2(c) and Fig. 3(c) show that 
the two additional similar objects enter from the x axis, 
forming a rhombus configuration and A±2 increase fur- 
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ther. Note, however, that the exact four-fold symmetry 
is not attained at this rotation frequency. 

3.2 Polar phase 

We now go on discussing the polar (PO) phase charac- 
terized by I /I =0 and \Q\ ^ 0. We take up the parameter 
values /3i = — 1 and /33 = 1 as a representative point for 
PO as shown in (c) in Fig. 1. Figure 5(a) displays the 
PO phase where \Aq\ ^ and all the others vanish, cor- 
respondingly 1/1=0 and |0| = 1 throughout the system. 
This is genuine PO configuration. 

Under rotation at = 0.18 shown in Fig. 5(b) the two 
HQVs emerge from the outside along the x axis. Namely 
while 1 remains intact, A±2 components mainly carry 
two HQVs. It is seen from Fig. 5(b) that |A_2| exhibits 
the suppressions on the x axis, which are compensated by 
A2 component. The OP can be described locally around 
one of those singular points by the half quantum vor- 
tex form: A(k) = exp(i(/7/2) exp(-i((^/4)F^)[y2,2(k) + 
i2,-2(k)] = y2,2(k) + e*'^r2,-2(k). Correspondingly at 
the core regions the ferromagnetic component is induced 
(see I/I and |6| in Fig. 5(b)). It should be noted here that 
since there is non-vanishing component of at those 
cores, those vortices are not true half-quantum vortex. 

At O = 0.25 it is seen from Fig. 5(c) those HQV-like 
objects increase in number which is now four. Moreover a 
new feature at this higher rotation emerges: four singular 
points in Aq are seen in Fig. 5(c). Thus four HQV-like 
objects and four singular cores in Aq coexist in the sys- 
tem, each of which are paired up as seen from |/| and |6| 
in Fig. 5(c). Note that those two kinds of objects show 
only non- vanishing /^ component. There is no fx and 
fy components throughout the system because A±i are 
completely absent. 

3.3 Cyclic phase 

We take up the values appropriate for CY; /3i = 
1, /?3 = —1 (sec (d) in Fig. 1). The following re- 
sults do not depend on these values and are generic in 
the cyclic region. By solving the GL equations, we ob- 
tain the stable CY in the presence of a trap. Among 
various CY forms derived from the canonical CY de- 
scribed by (i, 0, a/2, 0, Vy^ , we stabilize the particular CY 
phase, called CY-z: (1, 0, 0, \/2, 0)"^, which is obtained 
by exp(i cos-i(l/V3)J^y) exp(-i(7r/4)j;)(i, 0, v^, 0, if 
through rotations of the canonical CY. This CY-z will 
be seen to support the 1/3-vortex later. As seen from 
Fig. 6(a) CY-z is dominant in the central region because 
I/I = and simultaneously |0| = 0. The surface region 
consists of FM or PO phases, which are intermingled. 
Note that FM is advantageous in the surrounding re- 
gion because as mentioned the system can gain in the 
condensation energy by tuning the direction of /. This 
OP texture differs completely from the spinor F=2 BEC 
where the CY-z extends all the way to the surface of the 
cloud^^) because of different boundary conditions. We 
emphasize that wherever OP varies spatially, the gra- 
dient coupling terms fgrad in eq. (4) inevitably induce 
other components. 

The gap structure of CY consists of 8 point nodes; In 



canonical CY: (i, 0, v^, 0, z)"^ the point nodes are situ- 
ated at the 8 directions (±1, ±1, ±1) given by 8 corners 
of the inscribed cube inside the Fermi sphere. In CY-z 
the inscribed cube is rotated so that the (1,-1,1) di- 
rection lies now parallel to the z axis. Thus this CY-z 
is three-fold symmetric with respect to combined gauge 
transformation and rotations about the z axis. This is 
the origin for the possible existence of a 1/3-vortex in 
CY-z as discussed below. The mass current given by cq. 
(9) flows spontaneously circularly around the center of 
the trap and gradually decreases farther away from the 
origin. 

By rotating CY-z we can create the 1 /3- vortex whose 
main structure is described by (e*'^, 0, 0, v^, 0)-^ in the 
central region as seen from Fig. 6(b). Note that at the 
center the A^ component vanishes whereas the A-\ com- 
ponent is non-zero. Simultaneously the other components 
are induced around the surface region. Combination of 
the three-fold symmetry around the 2;-axis for CY-z with 
an additional gauge transformation leads to the 1/3- 
vortex form. Namely, CY-z: y2,2(k) + V2F2 _i(k) trans- 
forms into F2,2(k)e2'"^-|-\/2y2,-i(k)e-*'^ under a rotation 
of an angle (j) around the z-axis, which can be rewritten 
as e-*'^'(Y2,2(k)e3*'^ -h V2r2,_i(k)) = e-*'^/3(y2,2(k)e*'^ + 
\/2F2,-i(k)) by identifying = (/j/3. After the gauge 
transformation by (p/3 this becomes the 1/3-vortex form 
(e*''', 0, 0, v^, 0)"^. The vortex center is dominated by 
the A-i component where the A2 component vanishes 
due to the phase singularity at the core as seen from 
Fig. 6(b). Thus the core of the 1/3-vortex is ferromag- 
netic. The nodal structure of the ferromagnetic core re- 
gion is characterized by a line node on the equator of the 
Fermi sphere in addition to point nodes at both poles 
because the corresponding basis function is y2,-i(k) oc 

We also calculated the mass current for the 1 / 3- vortex 
and the line integral of the velocity field § v^dLp along 
various closed paths around the center. We find that the 
line integral yields approximately (l/3)(/i/2m) when cir- 
cling near the center, implying 1/3 quantization. How- 
ever, this conclusion is not exact because the system 
is inhomogeneous, a situation different from the infinite 
system where the quantization is exact if the line integral 
is taken along a path far away from the vortex line. 

It is remarkable that this 1/3-vortex is the most sta- 
ble vortex among various possible ones which we have 
solved in order to compare their energies. The critical ro- 
tation frequency is found to be = 0.21w at T = O.STg 
beyond which a single 1/3-vortex becomes energetically 
stable compared to the CY-z. 

4. Conclusions and summary 

We have examined theoretically textures and vortices 
for d-wave superfluids produced by a Feshbach resonance 
of neutral Fermionic atomic gases. Our arguments are 
based on Ginzburg-Landau framework which is general 
enough to capture the basic and fundamental properties 
of rf-wave superfluids confined by a harmonic trap be- 
cause GL functional is determined by symmetry based 
arguments; C/(l) gauge symmetry and 50(3) spatial ro- 
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Fig. 5. (Color online) Contour maps of \A2\ ~ |A_2|, \f\, and |e| in PO phase. /3i = -1 and /33 = 1. (a) H = 0, (b) H = 0.18, and (c) 
n = 0.25. 
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Fig. 6. (Color online) Contour maps of IA2I ~ |A_2|, |/|, and |e| in CY phase. /3i = 1 and ^3 = -1. (a) C = and (b) C = 0.22. 



tat ion group. 

We exhausts all three possible stable phases; ferromag- 
netic, polar and cyclic phases, and show the possible tex- 
tures and vortices for each phase at rest and under rota- 
tion. We hope that this study helps to identify the d-wave 
superfluidity experimentally. 

In particular, the 1/3- vortex characteristic to CY is 
intriguing because this is energetically stable over other 
possible vortices. A possible way to realize it experimen- 
tally as follows. The initial preparation of CY-z itself is 
not difficult where the appropriate ratio (1 : \/2) in IA2I 
and \A-i \ is needed. Since these two components are dif- 
ferent orbital angular momentum states, the population 
in the components can be adjusted using Raman tran- 
sitions, similarly as in F=2 spinor BEC^^'^^^ or by us- 
ing RF transition in the presence of Zeeman splitting of 
the 5 components under an applied field as was done by 
Hirano's group. Furthermore, once CY-z is prepared, 
it is possible to imprint the unit phase winding in the 
A2 component using a circular polarized Raman transi- 
tion. It may be difficult to use the optical spoon 
method to create the vortex, which was utilized to pro- 
duce vortices in scalar BECs.^^-' 
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Appendix 



+ i^i[(|aiB22p + 1925331' + lasBiiP + |aii333|2 + ISa^nl^ + 193^32 



In general, the contractions of the tensor {diB jk)* {d\B f^u) with the rank six are given by the following three forms: 
{diBjk)*{diBjk), {diBjk)*{dkBij), and (djBij)* (dkBik). The linear combination of these invariants gives rise to the 
gradient term of the GL functional: 

fgrad = Ki{diBjkT {diB^k) + K^idiBjk)* {dkBij) + K^{djB,,Y {dkBik) 

= {K, +K2 + K3) [|9iBn|' + 1925221' + 1535331'] 

'' ' " '' ' ' ^25i2|' + |93i>23| 

522 

+ 2(|9i523|' + |a253i|2 + |a35i2p)] 
+ /f2[(9l523)*(a253l) + (9253l)*(a35i2) + (935l2)* (91523) + (9l 533 )* (93531 ) + (^2511 )* (9i5i2) 

+ (^3522)* (92523) + (935n)*(9i53l) + (9i522)*(925i2) + (92533)* (93523) + C.C.] 
+ i^3[(a2523)*(9i53l) + (9353l)*(925i2) + (9i5i2)* (93523) + (9i5ii)*(9353l) + (92522)*(9i5i2) 

+ (93533)* (92523) + (93533)* (9i 531) + (9i5n)*(925i2) + (92 522)* (93523) + C.C.] . (14) 

On the other hand, the GL functional can be also written in terms of the OP Ak with k = —2,-1,0,1,2. In 
particular, the gradient term is rewritten as 

f, = -iV^C2(/i^;F)'i(n:^(k)>l,,(k)r2:,(k)y2,i(k))k(a^Afe)*(9,A0 (15) 
where Np is the density of states at the Fermi level. The coefficient C2 is given by 

where a;„ is the Matsubara frequency. Therefore 

Near T — > Tc, the above GL functional becomes 

farad = i^o(>^M(k)>l,,(k)r2%(k)y2,i(k))fc(a.Afc)*(9,^,) (18) 

where we define Kq = 7({3)Np{hvF)'^/4:8Tr'^kgT^. Among various combinations of the above indices i,j,k,l, the 
non-zero contributions are classified into three cases (A), (B), and (C): 
(A) i = j,k = l 

fgrad =^^0(1^1. I' |l'2,fe |') l^, Afe P 

=Ko[{\Y,,o\^\Y^,o\')\doAo\' + (|Fi,o|'|i1^,il') (l^oAip + 190^-1^ 
+ {\Yuo\''\Y2,2\^) {\d0A2f + |9o^-2|') + {\Yi,if\Y2,of) {\d+Ao\^ + |9_Ao|') 
+ (|yi,ir|y2,i|') (|S+^i|' + \d+A_^f + |9_Ai|2 + |9_A_i|2) 

+ (|n,i|'|l2,2|') ild+A^f + |9+A_2|' + |9_A2|' + \d-A_2f)] (19) 

where (|Yi,o|'|>2,o|') = 11/7, etc. The other cases (B) i = j ± 1, k = I T 1 and (C) i = j ± 2, A; = Z =F 2 are treated 
similarly. By summing up those terms (A), (B), and (C) we obtain 

fgrad ~ [ll|So^o|' + 9 {\daA^\^ + \doA_^\^) + 3 (|9oA2|2 + |9oA_2H + 5 (|9+Ao|' + 1^-^01') 

+ 6 (|9+Ai|2 + |9+A_i|2 + |9_Ai|2 + \d-A-^\^) + 9 (19+^21" + |9+A_2|' + |9_A2|' + |9_A_2|') 
+ {^^^[{^oA^f{^^Ao) + (9_^o)*(9o^-i) - {dM*{d+A_,) - (9_^i)*(9oAo)] 

+ 3y2[(9oA2)*(9+Ai) + (9_A_i)*(9oA_2) - (9oA_i)*(9+A_2) - (9_A2)*(9oAi)] 

+ 6(9_Ai)*(9+A_i) + 2V6[(9_A2)*(9+Ao) + (9_Ao)*(9+A_2)] + C.C.}] . (20) 
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This is further rewritten as 

ll|5,^oP + 9 {\d,A,\^ + + 3 {\d,A2\^ + \d,A_2\^) + 5 {\d,Ao\^ + \dyAo\^) 



' 7 



+ 6 {\d^Aif + \dyAi\'' + + \dyA_,\'') + 9 {\d^A2f + \dyA2\'' + \d,A_2f + \dyA_2\'') 

+ |^Pz^o)*{(5, - tdy)A_,} + {{d, + zdy)Aor{d,A_,) 

- {d,A,y{{d, - idy)Ao} - m + idy)Air{d,Ao)] 
+ 3p,^_i)*{(a, - idy)A-2} + {(9x + idy)A-iy{d,A-2) 

- {d,A2r{{dx - idy)A,} - {{dx + idy)A2}*{d,Al) - {{dx + idy)A,y{{dx - idy)A-,}] 



- V^[{(Sx + idy)A2}*{id^ - idy)Ao} + {{d^ + idy)Ao}*{{d^ - idy)A_2}] + C.C. 
where do = and d± = T{dx T idy)/y/2. After rearranging it, we finally obtain 



(21) 



fgrad =^|ll 1 91^32.2 .^2 |^ + |(9yA3j,2_r2 |^ + \dzA:iz2_r2\' 



+ 9 [\d,A,y\^ + \dyAy,\^ + |5,A,,p + \dM^ + \dv^.v? + \d.^yz?\ 

l^y^zxl^ + |f^2^Xl/P + |f^X^3/2_22| + \dyAz2_^2\ + 1 9^ j4j;2 _y2 j^] 
+ WxAy,r{dyA,^) + {dyA,^y{d,A^y) + (9, ) ^ (9, Aj,, ) 

+ {dyAy,y{d^A,^) + {d,A,^r{dyA^y) + {d ^A^yY {3 ,Ay,) + C.C] 
- 2\/3[(axA3y2_^2)*(9^A^^) + {dyA^^2_^2)*(dxA^y) + (9^ ^3^2 ) * (9y Ay ^ ) 

+ {d,A^y2_,2)*{d^A,^) + {d^A^,2_,2)*(dyA^y) + (dyA^^2_,2y{d,Ay,) + C.C.]} (22) 

where we transform the basis of A from the z component of the orbital angular momentum to the direction of the 
rectangular coordinates. 

In order to confirm that eq. (22) is equivalent to the gradient form eq. (14) written in terms of we first note 
that the transformation between B and A is given as 



-B33, A'ix2_j.2 



A^2_y2 — y — (iJll — B22) 



A:r.,i — \ — -Bio, A 



g-Bll, Asy2_r2 



Ay2_z2 



15 



12) ^yz 



15 



^^22, 



— {B22 - B33) , A,2_,2 = W — (S33 - Sn) 



15 ^ 



(23) 



By substituting B in place of A, we obtain 

2 
7 
12 
' 35 



fgrad 



Ko [\diBi,f + \d2B22f + \d3B3sf] 

Ko [\diB,2\^ + 1925231' + 1935311' + \dlB3l\^ + \d2Br2? + 1935231'] 



+ ^i^o[(|9i522r + 1925331' + |935n|' + |9i533|' + \d2Bn\^ + |93522|') 

+ 2 (|9l523|' + 192531 1' + |935i2|')] 
+ ^i^o[(5l523)*(9253i) + (9253i)*(935i2) + (935i2)*(9i523) 

+ (9i533)*(9353i) + (925n)*(9i5i2) + (SgSas)* (93523) 

+ (935n)*(9i53i) + (9i522)*(925i2) + (92533)* (93523) + C.C] 

+ ^ifo[(92523)*(9i53l) + (9353l)*(925i2) + (9i5i2)*(93523) 
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+ id3B33)*{dlB3l) + {dlBii)*{d2Bi2) + (52B22)* (93523) + C.C.] (24) 

where we use the traceless property of B. This is indeed the gradient GL functional eq. (14) written in terms of B 
above. The coefficients are 



2Ki =K2 = K3 = {4/35)Ko. 



(25) 
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